# Code from ~/Final_Code_Replication_Final/01_Appendix/Appendix K/

# Code for Table K1: Biprobit

rm(list = ls())

library('AER')
library('stargazer')
library('ivpack')
library('dplyr')
library('GJRM')
library('data.table')

## ---------------------------------------
## Load Data and Functions
## ---------------------------------------
load("../data0812.RData")

directory <- "../functions/"
functions <- list.files(directory)  
loadfunctions <- sapply(functions, FUN = function(x)source(paste0(directory, x)))


## ---------------------------------------
# Construct Instrument
## ---------------------------------------
data0812 <- constructIV(data0812)

## ---------------------------------------
#  Select Last Case Before Election
## ---------------------------------------
data0812 <- lastCase(data0812)

## ---------------------------------------
# convert probabilistic to binary voter file measures
## ---------------------------------------

data0812$vote2012 <- as.numeric(data0812$vote2012 > 0.80)
data0812$vote2008 <- as.numeric(data0812$vote2008 > 0.80)
data0812$regis_before <- as.numeric(data0812$regis_before > 0.80)

## ---------------------------------------
## Covariates
## ---------------------------------------

out.eq <- vote2012 ~ pti + as.factor(court_time1) + as.factor(court_time2) + 
  as.factor(court_dow) + as.factor(court_shift) + as.factor(totOGS2) + 
  age_2012 + I(age_2012^2) + Female + as.factor(race) + vote2008 + 
  as.factor(noteli08) + regis_before + as.factor(any_drug_2) + 
  as.factor(any_violent_2) + as.factor(fire_arms_2) + as.factor(any_rob_2) + 
  as.factor(any_dui_2) + as.factor(prior_offender_2)

treat.eq <- pti ~ judgeiv + 
  as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + 
  as.factor(court_shift) + as.factor(totOGS2) + age_2012 + 
  I(age_2012^2) + Female + as.factor(race) + vote2008 + as.factor(noteli08) + 
  regis_before + as.factor(any_drug_2) + as.factor(any_violent_2) + 
  as.factor(fire_arms_2) + as.factor(any_rob_2) + as.factor(any_dui_2) + 
  as.factor(prior_offender_2)

f.list <- list(treat.eq, out.eq)
mr <- c("probit", "probit")
bvp <- gjrm(f.list, data = data0812, Model ="B", margins = mr)

data0812_1 <- data0812_2 <- data0812
data0812_1$pti <- 1
data0812_2$pti <- 0

out1 <- predict(bvp, eq = 2, newdata = data0812_1, type = 'response')
out0 <- predict(bvp, eq = 2, newdata = data0812_2, type = 'response')
outF <- mean(out1) - mean(out0)
outF

## ---------------------------------------
## Standard Error: Bootstrap
## ---------------------------------------
sims <- 500

library(parallel)
library(doParallel) #contains foreach
cl <- detectCores()/2

registerDoParallel(cl) #create backend

out <- foreach(j = 1:sims) %dopar% {
  set.seed(j)
  index <- sample(1:nrow(data0812), nrow(data0812), replace = TRUE)
  
  out.eq <- vote2012 ~ pti + as.factor(court_time1) + as.factor(court_time2) + 
    as.factor(court_dow) + as.factor(court_shift) + as.factor(totOGS2) + 
    age_2012 + I(age_2012^2) + Female + as.factor(race) + vote2008 + 
    as.factor(noteli08) + regis_before + as.factor(any_drug_2) + 
    as.factor(any_violent_2) + as.factor(fire_arms_2) + as.factor(any_rob_2) + 
    as.factor(any_dui_2) + as.factor(prior_offender_2)
  treat.eq <- pti ~ judgeiv + 
    as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + 
    as.factor(court_shift) + as.factor(totOGS2) + age_2012 + 
    I(age_2012^2) + Female + as.factor(race) + vote2008 + as.factor(noteli08) + 
    regis_before + as.factor(any_drug_2) + as.factor(any_violent_2) + 
    as.factor(fire_arms_2) + as.factor(any_rob_2) + as.factor(any_dui_2) + 
    as.factor(prior_offender_2)
  
  f.list <- list(treat.eq, out.eq)
  mr <- c("probit", "probit")
  bvp <- gjrm(f.list, data = data0812[index, ], Model="B", margins = mr)
  
  data0812_1 <- data0812_2 <- data0812[index, ]
  data0812_1$pti <- 1
  data0812_2$pti <- 0
  
  out1 <- predict(bvp, eq = 2, newdata = data0812_1, type = 'response')
  out0 <- predict(bvp, eq = 2, newdata = data0812_2, type = 'response')
  mean(out1) - mean(out0)
}

stopImplicitCluster()

cat("\nPrinting Table K1: Main Effects (Bivariate Probit)...\n")
print(paste("Estimate: ", round(outF, 3), "   Std. Error: ", round(sd(do.call('rbind', out)), 3)))
      